#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>

#define MAX_POLY_LEN 100    // 最大多项式系数长度
#define MAX_STR_LEN 2048    // 最大字符串长度

// ========================= 霍纳多项式结构体（复刻HornerPolynomial） =========================
typedef struct {
    double* coefficients;  // 多项式系数（高次到低次，如ax³+bx²+cx+d → [a,b,c,d]）
    int degree;            // 多项式次数（系数长度-1）
} HornerPolynomial;

// 创建霍纳多项式
HornerPolynomial* horner_poly_create(double* coeffs, int len) {
    HornerPolynomial* poly = (HornerPolynomial*)malloc(sizeof(HornerPolynomial));
    poly->degree = len - 1;
    poly->coefficients = (double*)malloc(len * sizeof(double));
    memcpy(poly->coefficients, coeffs, len * sizeof(double));
    return poly;
}

// 释放霍纳多项式内存
void horner_poly_destroy(HornerPolynomial* poly) {
    if (!poly) return;
    free(poly->coefficients);
    free(poly);
}

// 霍纳法则求值（核心）
double horner_poly_evaluate(const HornerPolynomial* poly, double x) {
    double result = 0.0;
    for (int i = 0; i <= poly->degree; i++) {
        result = result * x + poly->coefficients[i];
    }
    return result;
}

// 霍纳多项式转字符串（复刻__str__）
char* horner_poly_to_str(const HornerPolynomial* poly) {
    char* str = (char*)malloc(MAX_STR_LEN * sizeof(char));
    str[0] = '\0';
    bool first = true;
    int n = poly->degree + 1;

    for (int i = 0; i <= poly->degree; i++) {
        double coeff = poly->coefficients[i];
        if (fabs(coeff) < 1e-10) continue;  // 系数为0跳过

        // 符号处理
        if (!first) {
            strcat(str, "+");
        }
        first = false;

        // 系数部分
        if (fabs(coeff) != 1.0 || (n - i - 1) == 0) {
            char coeff_str[64];
            sprintf(coeff_str, "%.10f", coeff);
            strcat(str, coeff_str);
            strcat(str, "*");
        }

        // 次数部分
        int order = n - i - 1;
        if (order > 0) {
            strcat(str, "x");
            if (order > 1) {
                char order_str[16];
                sprintf(order_str, "**%d", order);
                strcat(str, order_str);
            }
        } else {
            strcat(str, "1");
        }
    }

    // 空多项式返回"0"
    if (first) {
        strcpy(str, "0");
    }
    return str;
}

// ========================= 分段多项式结构体（复刻PeacewisePolynomials） =========================
typedef struct {
    double* konts;                // 分段节点（x坐标）
    int kont_count;               // 节点数量
    HornerPolynomial** polynomials; // 分段多项式数组
    int poly_count;               // 多项式数量
} PeacewisePolynomials;

// 创建分段多项式
PeacewisePolynomials* peacewise_poly_create(double* konts, int kont_cnt, HornerPolynomial** polys, int poly_cnt) {
    PeacewisePolynomials* pwp = (PeacewisePolynomials*)malloc(sizeof(PeacewisePolynomials));
    pwp->konts = (double*)malloc(kont_cnt * sizeof(double));
    memcpy(pwp->konts, konts, kont_cnt * sizeof(double));
    pwp->kont_count = kont_cnt;
    pwp->polynomials = polys;
    pwp->poly_count = poly_cnt;
    return pwp;
}

// 释放分段多项式内存
void peacewise_poly_destroy(PeacewisePolynomials* pwp) {
    if (!pwp) return;
    free(pwp->konts);
    for (int i = 0; i < pwp->poly_count; i++) {
        horner_poly_destroy(pwp->polynomials[i]);
    }
    free(pwp->polynomials);
    free(pwp);
}

// 分段多项式求值（复刻__call__）
double peacewise_poly_evaluate(const PeacewisePolynomials* pwp, double x) {
    for (int i = 0; i < pwp->poly_count; i++) {
        double x0 = pwp->konts[i];
        double x1 = pwp->konts[i + 1];
        if (x >= x0 - 1e-10 && x <= x1 + 1e-10) {  // 浮点精度容错
            return horner_poly_evaluate(pwp->polynomials[i], x - x0);
        }
    }
    return NAN;  // 超出区间返回NaN
}

// 分段多项式转字符串（复刻__str__）
char* peacewise_poly_to_str(const PeacewisePolynomials* pwp) {
    char* str = (char*)malloc(MAX_STR_LEN * pwp->poly_count);
    str[0] = '\0';

    for (int i = 0; i < pwp->poly_count; i++) {
        char* poly_str = horner_poly_to_str(pwp->polynomials[i]);
        char segment_str[MAX_STR_LEN];
        sprintf(segment_str, "%s,[%.0f,%.0f]\n", poly_str, pwp->konts[i], pwp->konts[i + 1]);
        strcat(str, segment_str);
        free(poly_str);
    }

    return str;
}

// ========================= 辅助函数 =========================
// 数据点排序（按x坐标升序，复刻sorted + getx）
void sort_values(double** values, int rows) {
    for (int i = 0; i < rows - 1; i++) {
        for (int j = 0; j < rows - i - 1; j++) {
            if (values[j][0] > values[j + 1][0]) {
                // 交换两行
                double* temp = values[j];
                values[j] = values[j + 1];
                values[j + 1] = temp;
            }
        }
    }
}

// ========================= 三次样条插值核心算法（μ数列法） =========================
PeacewisePolynomials* cubic_spline_interpolate(double** values, int rows) {
    if (rows < 2) {
        fprintf(stderr, "数据点数量至少为2\n");
        return NULL;
    }

    // 步骤1：排序数据点
    sort_values(values, rows);

    // 步骤2：提取x、y数组，计算h数组
    int n = rows - 1;
    double* x = (double*)malloc(rows * sizeof(double));
    double* y = (double*)malloc(rows * sizeof(double));
    double* h = (double*)malloc(n * sizeof(double));

    for (int i = 0; i < rows; i++) {
        x[i] = values[i][0];
        y[i] = values[i][1];
    }
    for (int i = 0; i < n; i++) {
        h[i] = x[i + 1] - x[i];
    }

    // 步骤3：求解μ和z数组
    double* mu = (double*)calloc(n, sizeof(double));
    double* z = (double*)calloc(rows, sizeof(double));

    for (int i = 1; i < n; i++) {
        double g = 2 * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];
        mu[i] = h[i] / g;

        double temp = 3 * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1]) + y[i - 1] * h[i]);
        z[i] = (temp / (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g;
    }

    // 步骤4：求解系数a、b、c
    double* a = (double*)malloc(n * sizeof(double));
    double* b = (double*)calloc(rows, sizeof(double));
    double* c = (double*)malloc(n * sizeof(double));

    for (int j = n - 1; j >= 0; j--) {
        b[j] = z[j] - mu[j] * b[j + 1];
        c[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (b[j + 1] + 2 * b[j]) / 3;
        a[j] = (b[j + 1] - b[j]) / (3 * h[j]);
    }

    // 步骤5：构建分段霍纳多项式
    HornerPolynomial** polys = (HornerPolynomial**)malloc(n * sizeof(HornerPolynomial*));
    for (int i = 0; i < n; i++) {
        // 三次多项式系数：[a[i], b[i], c[i], y[i]]（对应ax³+bx²+cx+d）
        double coeffs[4] = {a[i], b[i], c[i], y[i]};
        polys[i] = horner_poly_create(coeffs, 4);
    }

    // 步骤6：创建分段多项式
    PeacewisePolynomials* result = peacewise_poly_create(x, rows, polys, n);

    // 释放临时内存
    free(x);
    free(y);
    free(h);
    free(mu);
    free(z);
    free(a);
    free(b);
    free(c);

    return result;
}

// ========================= 测试代码 =========================
int main() {
    // 测试数据：[[1,1],[2,2],[3,1],[4,-2],[5,1],[6,2],[7,1],[8,-2],[9,1]]
    int rows = 9;
    double** values = (double**)malloc(rows * sizeof(double*));
    values[0] = (double[]){1.0, 1.0};
    values[1] = (double[]){2.0, 2.0};
    values[2] = (double[]){3.0, 1.0};
    values[3] = (double[]){4.0, -2.0};
    values[4] = (double[]){5.0, 1.0};
    values[5] = (double[]){6.0, 2.0};
    values[6] = (double[]){7.0, 1.0};
    values[7] = (double[]){8.0, -2.0};
    values[8] = (double[]){9.0, 1.0};

    // 三次样条插值计算
    PeacewisePolynomials* cubic = cubic_spline_interpolate(values, rows);
    if (!cubic) {
        fprintf(stderr, "插值失败\n");
        return -1;
    }

    // 测试求值
    double x1 = 1.5;
    printf("cubic(%.1f) = %.10f\n", x1, peacewise_poly_evaluate(cubic, x1));

    // 打印分段多项式
    char* cubic_str = peacewise_poly_to_str(cubic);
    printf("%s", cubic_str);

    // 释放内存
    free(cubic_str);
    peacewise_poly_destroy(cubic);
    free(values);

    return 0;
}